/*
 * phase_lattice_mkl.hpp
 *
 * Copyright 2009-2012 Karsten Ahnert
 * Copyright 2009-2012 Mario Mulansky
 *
 * Distributed under the Boost Software License, Version 1.0.
 * (See accompanying file LICENSE_1_0.txt or
 * copy at http://www.boost.org/LICENSE_1_0.txt)
 */

#ifndef PHASE_LATTICE_MKL_HPP_
#define PHASE_LATTICE_MKL_HPP_

#include <cmath>

#include <mkl_blas.h>
#include <mkl_vml_functions.h>
#include <boost/array.hpp>

template <size_t N>
struct phase_lattice_mkl {
  typedef double value_type;
  typedef boost::array<value_type, N> state_type;

  value_type m_epsilon;
  state_type m_omega;
  state_type m_tmp;

  phase_lattice_mkl()
    : m_epsilon(6.0 / (N * N))  // should be < 8/N^2 to see phase locking
  {
    for (size_t i = 1; i < N - 1; ++i)
      m_omega[i] = m_epsilon * (N - i);
  }

  void inline operator()(const state_type& x, state_type& dxdt, const double t) {
    const int n = x.size();

    dxdt[0] = m_omega[0] + sin(x[1] - x[0]);

    vdSub(n - 1, &(x[1]), &(x[0]), &(m_tmp[0]));
    vdSin(n - 1, &(m_tmp[0]), &(m_tmp[0]));
    vdAdd(n - 2, &(m_tmp[0]), &(m_tmp[1]), &(dxdt[1]));
    vdAdd(n - 2, &(dxdt[1]), &(m_omega[1]), &(dxdt[1]));

    dxdt[N - 1] = m_omega[N - 1] + sin(x[N - 1] - x[N - 2]);
  }
};

#endif /* PHASE_LATTICE_MKL_HPP_ */
